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ABSTRACT 

We explore the possibility that the observed eccentricity distribution of extrasolar planets arose 
through planet-planet interactions, after the initial stage of planet formation was complete. Our 
results are based on ~ 3250 numerical integrations of ensembles of randomly constructed planetary 
systems, each lasting 100 Myr. We find that for a remarkably wide range of initial conditions the 
eccentricity distributions of dynamically active planetary systems relax towards a common final equi- 
librium distribution, well described by the fitting formula dn cx e exp[— i(e/0.3) 2 ]c?e. This distribution 
agrees well with the observed eccentricity distribution for e > 0.2, but predicts too few planets at 
lower eccentricities, even when we exclude planets subject to tidal circularization. These findings sug- 
gest that a period of large-scale dynamical instability has occurred in a significant fraction of newly 
formed planetary systems, lasting 1-2 orders of magnitude longer than the ~ 1 Myr interval in which 
gas-giant planets are assembled. This mechanism predicts no (or weak) correlations between semima- 
jor axis, eccentricity, inclination, and mass in dynamically relaxed planetary systems. An additional 
observational consequence of dynamical relaxation is a significant population of planets (> 10%) that 
are highly inclined (> 25°) with respect to the initial symmetry plane of the protoplanetary disk; this 
population may be detectable in transiting planets through the Rossitcr-McLaughlin effect. 
Subject headings: planetary systems - planetary systems: formation - planets and satellites: general 



1. INTRODUCTION 

Our understanding of the origins of planetary systems 
has been revolutionized in the last decade by the discov- 
ery of over 250 extrasolar planets. These discoveries have 
vastly broadened our appreciation of the diversity of pos- 
sible planetary systems, and raise a number of challenges 
for theories of planet formation. 

One of the most important of these is the problem 
of large eccentricities. The median eccentricity of the 
known extrasolar planets is ~ 0.2 (even before discard- 
ing those planets at small semimajor axes whose orbits 
have likely been circularized by tidal forces), a factor of 
two larger than that of any solar-system planet except 
Mercur y. The largest kno wn eccentricity is 0.93, for HD 
80606b (jNaef et al.ll2001h . The large eccentricities are a 
major concern because the nearly circular, coplanar or- 
bits of solar-system planets have been one of the strong 
arguments that the planets formed from a disk since the 
time of Kant and Laplace. The mechanism by which the 
extrasolar planets acquired their eccentricities, and why 
the eccentricities of the planets in the solar system are so 
much smaller than those of the known extrasolar planets 
(the "eccentricity problem" ) are still unknown. The res- 
olution of the eccentricity problem, and the wider ques- 
tion of understanding the distributions and correlations 
of the dynamical and physical parameters of planets, are 
key milestones on the road towards a comprehensive the- 
ory of planet formation. 

The process of planetary system formation is notori- 
ously difficult to study theoretically. An oversimplified 
but useful approach is to divide it into two stages, based 
on the importance of the effects of gas and the protoplan- 
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etary disk on the growing planets. Stage 1 lasts a few 
Myr, until the dissipation of the gaseous protoplanetary 
disk. Irrespective of the exact formation mechanism, at 
the end of this stage the gas giant planets should have 
accreted their massive envelopes — they must do so be- 
fore the gas disk disappears — and some (very uncer- 
tain) number of smaller solid planctcsimals, planetary 
embryos, and planets should also be present. 

Stage 2 lasts from the end of Stage 1 until the present. 
In Stage 2 the evolution is driven primarily by gravi- 
tational interactions and collisions between the planets. 
These processes can lead to the ejection of planets into 
interstellar space or the Oort cloud, the consumption of 
planets by the host star, and collisions and mergers of 
planets. All of these reduce the number N p i of surviving 
planets; presumably as N p i declines, the system gradu- 
ally evolves to a more and more stable state, and the 
timescale for future evolution is always comparable to 
the present age. While there is a vast literature on the 
origins of solar systems that focuses on Stage 1, investi- 
gations of Stage 2 evolution are surprisingly scarce. Yet 
Stage 2 is actually easier to explore than Stage 1 in most 
respects, since the processes in Stage 2 involve only the 
simplest possible physics: Newton's law of gravity and 
Newton's laws of motion. 

The first numerical explorati ons of Stage 2 wer e 
motivated by the suggest i on of iRasio fc Fordl (1 19961) . 
IWeidenschilling fc Marzarl (|1996h . and lLin fc Idal (|1997l ) 
that gravitational interactions between planets could be 
responsible for the large eccentricities of the extrasolar 
planets 3 . These numerical experiments demonstrated 

3 An additional motivation was the suggestion that hot Jupiters 
(massive planets on circular orbits with semimajor axes < 0.03 AU) 
might be formed by tidal circularization of high ly eccentric orbit s 
with initial semimajor axes of a few AU (e.g. IFord etHl feOOl); 
but, as the sample of extrasolar planets grew, this mechanism as 
originally proposed proved unable to explain the frequency of hot 
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that planet-planet scattering could excite eccentricities 
to the required levels. However, the simulations were 
somewhat unrealistic in that they only lasted ~ 10 4 - 
10 6 yr, shorter than the likely duration of Stage 1, and 
typically included only a few planets (usually two or 
three). Followup numerical investigations of the evo- 
lution of unstab le systems conta ining two identical gi- 
ant planets fe.g.. lFord et al.ll2001[ ) found them incapable 
of quantitatively matching the observed eccentricity dis- 
tribution, instead producing eccentricities smaller than 
th ose typically observed ( see also the analytic arguments 
in IGoldreich fc Sarill2003h . However, simulations of un- 
stable planeta ry systems conta ining two planets with un- 
equal masses dFord et al . 2003), or three or m ore pla nets 
(|Marzari fc Weidenschillingj l2002t iFord et al.ll2003h . do 
appear to be capable of exciting a sufficient fraction of 
planets to high eccentricities. 

Other works have examined sys tems with mor e plan - 
ets and over longer timescales. iLevison et al.l (|1998[ ) 
conducted a set of 16 long-term (1 to 16 Gyr) simula- 
tions of systems with ~ 10 planets, showing that sta- 
ble highly eccentric orbits are a possible outcome of 
long-t erm dynamical evolution. iPapaloizou fc Terqueml 
(|2001f ) reached a similar conclusion, although they only 
followed the evolution of their systems for < 10 4 or- 
bital periods, which is too short to distinguish Stage 
1 from Stage 2 evolution. They also argued, following 
iRasio fc Ford! (|1996l ). that dynamical evolution almost 
inevitably leads to ejections, which could contribute to a 
population of free-floating planets that may be detectable 
by larg e-scale microlensing s urvey s. In the largest study 
to date. lAdams fc Laugh lin (2003J) noted the importance 
of studying large ensembles of systems with similar ini- 
tial conditions, and observing the evolution of ensemble 
averages and distributions. They traced the evolution 
over 10 6 yr of four ensembles of 100 systems each, and 
the evolution over 10 7 yr of two ensembles of 50 systems, 
again concluding that dynamical mechanisms are capable 
of producing high eccentricities, as seen in the observa- 
tions. 

These studies suggest that Stage 2 evolution could play 
an important role in determining the numbers, masses, 
and orbital elements of extrasolar planets, but the short 
durations and small ensemble sizes of these integrations 
limit the impact of their conclusions. They also leave 
many questions unanswered; in particular, the relative 
importance of Stage 1 and Stage 2 in determining the 
properties of planetary systems is still largely unex- 
plored. An extreme (and unlikely) possibility is that 
the dynamical properties of planets are determined al- 
most entirely in Stage 2, so that planetary systems are 
in some kind of quasi-equilibrium state that is largely 
independent of the initial conditions set at the end of 
Stage 1 — just as the Maxwell-Boltzmann distribution 
of gas molecules in a room is determined entirely by the 
temperature and total gas mass, or the distribution of 
dark matter in a galactic halo appears to be almost in- 

Jupiters compared to planets a t larger distances. However, recent 
work by Naeasawa ct al. (2008), as well as the freque ncy of planet- 
star collisions observed in our simulations ('Section I2.3|l . indicate 
that planet scattering in systems with multiple planets, or in com- 
bination with the Kozai mechanism, might still plausibly explain 
the observed frequencies of hot Jupiters. This possibility deserves 
further investigation. 



dependent of how the halo formed (|Navarro "et~allfl99l . 
An equally extreme, and almost as unlikely, possibility is 
that the properties of planetary systems are determined 
entirely in Stage 1, so that the distribution of planetary 
systems looks the same at an age of a few Myr as it does 
after a few Gyr. 

This paper is the beginning of a systematic investiga- 
tion of the role of Stage 2 evolution in determining the 
properties of planetary systems, using accurate long-term 
(10 s yr or more) integrations of large (up to N sys = 1000) 
ensembles of mock planetary systems with up to 50 plan- 
ets per system. In terms of the natural metric for the 
computational effort required 

PP = (number of planets per system) x (1) 
(number of orbital periods) x 
(number of systems) 

the results presented here have PP = 5 x 10 12 , a fac- 
tor of 50 more than any previous exploration of Stage 2 
evolution. 

In this paper we focus on the dynamical evolution of 
the eccentricity distribution. In Sj2]we briefly describe the 
integrator, the selection of initial conditions, and simu- 
lation results. In $3] we discuss the distributions of ec- 
centricities of simulated systems, their classification as 
"active" , "inactive" and "partially active" , and a com- 
parison to observations. We quantify and justify the clas- 
sification criteria in SI SI discusses the inclinations of 
active systems, while $6] summarizes the results and dis- 
cusses their implications. 

2. SIMULATIONS 
2.1. The Integrator 

The principal challenge in following Stage 2 evolution 
is the computational intensity of the problem. Reliable 
studies require (i) accurate numerical integrations of N- 
body systems for at least ~ 10 s yr (our longer integra- 
tions show that relatively little evolution occurs between 
10 s and 10 9 yr); (ii) the ability to follow accurately high- 
eccentricity orbits and close encounters between planets 
(a challenge that is not present in long solar-system in- 
tegrations, where the planets are on well-separated, low- 
eccentricity orbits); (iii) large ensembles of planetary sys- 
tems, in order to make statistically significant predictions 
and explore the wide variety of possible outcomes. 

We overcome these challenges by the combination of 
an integration strategy tailored for the problem at hand, 
careful hand-optimization of key subroutines of the code, 
and the use of large computer clusters. To efficiently 
simulate thousands of planetary systems in a reasonable 
amount of time we have decided on a "mix-and-match" 
strategy, infrequently switching between the Bulirsch- 
Stoer and hybrid symplectic integration schemes de- 
scribed below, depending on the conditions present in 
the system being integrated. Our code makes use of the 
publicly availabl e integration routines of MERCURY6 4 
(|Cham bers 1999[), and extends them to support a high- 
level integration scheme as follows. 

Since the initial conditions of many of the simulated 
planetary systems produce dynamically active (numer- 
ous close encounters, frequent scatterings to high ec- 

4 Available at |http://www. arm. ac.uk/~jec /] 
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centricity orbits, etc.) and short-lived systems, the in- 
tegration of the first 10 6 yr is done using a high accu- 
racy (e = 10~ 12 ) conservative Bulirsch-St oer integrator 
C"BS2". [ChamberslH99llPress et al.lll992T) . 

For the subsequent much longer integration (10 6 - 
10 s yr), we switch to a hyb rid-symplectic ("HYBRID") 
integrator (|Chambersl lT999). This integrator accurately 
handles close encounters between planets, without the 
loss of symplectic properties. It accomplishes this by in- 
troducing a Hamiltonian splitting which, in the absence 
of close encounters, is equiva lent to the classical mixed - 
variable symplectic scheme (jWisdom fc Holmanl 1991). 
However, during a close encounter (defined as two plan- 
ets approaching to closer than a preset number of Hill 
radii, usually 3), the HYBRID algorithm, by means of 
a changeover function sensitive to mutual planet sepa- 
rations, redistributes the (now large) perturbation due 
to the encountering planets to the Keplerian part of 
the Hamiltonian. This keeps the remaining perturba- 
tion terms small, but makes the Keplerian part analyti- 
cally insoluble. The Keplerian part is therefore solved 
numerically to machine precision, using a high preci- 
sion Bulirsch-Stoer integrator. The complete solution 
is advanced in this manner to the end of the close en- 
counter, when the changeover function causes the HY- 
BRID scheme to again become equivalent to the compu- 
ta tionally efficient M VS. For details we refer the reader 
to lChamberd (fl999h . 

The HYBRID scheme handles close encounters with 
planets accurately, but can be susceptible to errors due 
to in adequate pericenter samp ling in highly eccentric or- 
bits (jRauch fc Holmanl Il999f ). To guard against this, 
during the symplectic integration phase the conservation 
of total energy and the individual orbital elements of 
planets in the system are continuously monitored. If 
the relative energy error averaged over the past 1000 
timesteps, of size h symp i, exceeds 10 , or if the dura- 
tion of pericenter passage (T ppn defined to be the time 
in which the true anomaly goes from — 7r/2 to +tt/2) of 
any planet become less than 5 timesteps, the code back- 
tracks At = 10000 timesteps and restarts the integration 
using the high-accuracy BS2 integrator (an "algorithm 
switch" ). The subsequent BS2 integration phase lasts for 
At = 30000h symp i, after which we return to the hybrid- 
symplectic integration with h symp i adjusted to l/15th of 
the smallest value of T pp in the system. We call this 
the "HYBRID/E" scheme since it is able to handle close 
pericenter passages and scatterings to high-eccentricity 
orbits, as well as close planetary encounters. 

The frequency of algorithm switches depends on the 
dynamical configuration of the system being integrated. 
As symplecticity is violated at each algorithm switch, it is 
desirable to keep these at a minimum. For the ensembles 
described in the next section, 60% of systems had no 
algorithm switches, 80% had fewer than 10, while only 
5% had more than 100 in 10 s yr of integration. 

We have tested the accuracy of this scheme by com- 
paring the results (energies, angular momenta, and time 
evolution of eccentricity) of three-body integrations of 
highly eccentric orbits using the BS2, HYBRID, and HY- 
BRID/E schemes (with more steps/orbit in the first two); 
typically, the errors of the HYBRID /E scheme are com- 
parable to those of HYBRID with ~ 2 times shorter 
timestep. We have also tested our algorithm by com- 



paring the results of 10 7 yr integrations of systems in 
the nlOslO ensemble (see Table [Hand the following sec- 
tion) with an integration done using the high-precision 
Bulirsch-Stoer scheme. The final distributions of eccen- 
tricities, inclinations, semimajor axes and masses were 
identical, within statistical error. 

2.2. Initial Conditions 

The selection of initial conditions would ideally be 
based on the predictions of Stage 1 planetary formation 
theory; unfortunately, the theory is still too crude to al- 
low this. We therefore picked the distributions of initial 
conditions (especially the eccentricities, see Figure [3]) for 
each ensemble in a largely arbitrary fashion, with only a 
minimal constraint that the planets begin in some sort 
of a disk. 

The ensemble definitions are detailed in Table 1. We 
chose semimajor axes uniformly in log(a), between a = 
0.1 and 100 AU. Similarly, we drew the masses from a 
distribution uniform in log(M), between M = 0.1 and 
10 Jupiter masses (th is is comparable to the obse rved 
distribution; e.g., see iTabachnik &: Tremaind 120021 and 
iMarcv et al.ll2005l ) . For all but one ensemble we drew the 
eccentricity e and the inclinations I from a Schwarzschild 
distribution 5 (jBinnev k. Tremaind [2008h : 

x f x^ \ 
dN = S(x; a x )dx = — ^ exp — — - ) dx (2) 

<?£ V 2(J xJ 

where x is either e or /, with o~ e and o~i as given in 
Table 1. If an eccentricity greater than 1 is drawn, the 
drawing is repeated until e < 1 is obtained. This effective 
truncation of the Schwarzschild distribution at e = 1 is of 
practical relevance in only one ensemble (cl0s40, having 
a e = 0.4), as S(x) falls off exponentially fast after the 
peak. Finally, the initial number of planets was either 3, 
10 or 50, depending on the ensemble. 

The planets were approximated by homogeneous 
spheres of density p = 1 gcm~ 3 . The central star was 
taken to have solar mass and radius. Planet-star colli- 
sions were allowed in all ensembles. For all but one en- 
semble (nlOslO) planet-planet collisions were allowed as 
well. Both planet-star and planet-planet collisions were 
assumed to be fully inelastic, resulting in momentum- 
conserving mergers with no fragmentation. We neglected 
all other effects, such as tidal dissipation and relativistic 
corrections. 

We have made no attempt to constrain our initial con- 
ditions to remove systems that are unstable over short 
timescales, e.g., by requiring that the planets be sepa- 
rated by some minimum number of Hill radii. The reason 
is that such systems are short-lived and hence consume 
a negligible fraction of our computing resources. 

For each ensemble, we constructed N sys realizations 
(planetary systems) which were then integrated for 
10 8 yr. This timespan corresponds to 3.2 • 10 9 and 10 5 
orbits at a = 0.1 and 100 AU, respectively. 

2.3. The average number of planets 

The time history of the average number of planets per 
system is plotted in Figure [T] The number of planets in 
all but one ensemble (c03s00, which contains only three 

5 Also known as the Rayleigh distribution. 
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TABLE 1 
Initial conditions 



Ensemble 


/(e) 


fW 


/(«) 








N pl 


Nsys 


Collisions 




Class 


c03s00 





S(7;0.05) 


U(loga; 


-1,2) 


U(logM 


-1,1) 


3 


500 


yes 


inactive 


nlOslO 


S(e;0.1) 


S(/;5.73) 


U(log a; 


-1,2) 


U(logM 


-1,1) 


10 


1000 


no 


P- 


active 


clOslO 


S(e;0.1) 


S(/; 5.73) 


U(log a; 


-1,2) 


TJflogM 


-1,1) 


10 


200 


yes 


P- 


active 


clOsOO 





S(/;0.05) 


U(loga; 


-1,2) 


U(logM 


-1,1) 


10 


150 


yes 


P- 


active 


cl0s30 


S(e;0.3) 


S(/;0.3) 


U(log a; 


-1,2) 


U(logM 


-1,1) 


10 


200 


yes 




active 


c50s05 


S(e;0.05) 


S(/;0.05) 


U(log a; 


-1,2) 


UYlogM 


-1,1) 


50 


200 


yes 




active 


cl0s40 


S(e;0.4) 


S(/;0.2) 


U(log a; 


-1,2) 


U(logAf 


-1,1) 


10 


500 


yes 




active 


c!0u80 


U(e; 0,0.8) 


S(/;3) 


U(log a; 


-1,2) 


U(logAf 


-1,1) 


10 


500 


yes 




active 



Note. — Definition of initial conditions for the ensembles of integrations. The columns (left to right) list 
the name of the ensemble, the distribution functions used for drawing the initial eccentricity /(e), inclination 
/(/) (degrees), semimajor axis /(a) (in AU), mass f(M) (in units of Jupiter mass), the initial number of 
planets N p i, and the number of realizations of the ensemble N sys . The column labeled Collisions specifies 
whether planet-planet collisions were allowed to occur during the simulation. The final column lists the 
ensemble classification according to the criteria of Section[3] All systems were integrated for 10 s yr. S(x; tr) is 
the Schwarzschild distribution (eq. [3J, V(x; x m i n , x max ) is the uniform distribution with x m i n < x < x max - 
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Fig. 1. — Average number of planets vs. time. Each curve shows 
the time evolution of (N p i) for systems of a particular ensemble, 
as marked by different symbols. 



planets) exhibits a rapid dropoff, starting on a dynamical 
(~-10 3 yr) and continuing on a secular timescale (10 3 - 
10 5 yr). A detailed breakdown, by mode of removal of 
planets from the system, is given in Table [2j In the 
first phase, the dominant mode of removal (the domi- 
nant "decay channel") is through planet-planet mergers 
as the randomly placed nearby planets, especially in the 
inner (a < 1 AU) regions, collide and merge. In the 
first 10 3 yr between 3% and 20% of all planets are lost 
to planet-planet collisions. After 10 3 yr there are only a 
few collisions, with none occurring in any of the ensem- 
bles after 10 6 yr. On t > 10 6 yr timescales, ejections to 
interstellar space become the dominant decay channel, 
with between 50% and 60% of planets being lost in this 



TABLE 2 

Number of planets lost to ejections, mergers, and 
stellar collisions 



Ensemble 


t - 


= 10^ 


yr 


t = 


= 10 b 


yr 


t = 


= 10* 


yr 




E 


M 


s 


E 


M 


s 


E 


M 


s 


c03s00 


0.0 


0.1 


0.0 


0.1 


0.1 


0.0 


0.1 


0.1 


0.0 


nlOslO 


0.5 


0.0 


0.3 


5.1 


0.0 


1.1 


6.2 


0.0 


1.3 


clOslO 


0.4 


0.5 


0.3 


1.0 


0.8 


1.6 


4.8 


0.8 


1.8 


clOsOO 


0.3 


1.2 


0.2 


3.0 


1.5 


1.3 


3.9 


1.5 


1.7 


cl0s30 


0.5 


0.7 


0.5 


4.2 


0.9 


1.9 


4.8 


0.9 


2.0 


cl0u80 


0.5 


0.5 


0.7 


1.2 


0.6 


2.3 


1.8 


0.6 


2.5 


cl0s40 


1.0 


0.7 


0.8 


4.3 


0.8 


2.1 


4.8 


0.8 


2.3 


c50s05 


5.7 


9.6 


5.1 


27.6 


9.7 


9.4 


28.9 


9.7 


9.5 



Note. — The average number of planets which by time t 
are lost to ejections (column E), mergers through planet-planet 
collisions (column M) and collisions with the star (column S). In 
ensemble nlOslO the values in column M are always 0, because 
planet-planet collisions are disallowed. 



way 6 . A further ~ 20% are lost to collisions with the 
star, usually as a result of gradual eccentricity excitation 
by a more massive planet. 

It is important to point out that our treatment of plan- 
ets on such collision orbits is not complete, as we neglect 
the dissipative tidal forces from the host that become 
significant when the pericenter is < 20 stellar radii. In 
particular, dissipative tides may circularize the planet 
orbit at small semimajor axis before it collides with the 
star. This mechanism may be responsible for some, but 
probably not most, of the hot Jupiters. Similarly, our 
treatment of ejected orbits is not complete, because we 
neglect the tidal forces from the Galaxy that become sig- 
nificant when the apocenter is > 10 3 AU. Galactic tides 
may cause some or even most planets on high-eccentricity 
orbits to end up in bound orbits of ~ 10 4 AU (analogous 
to the orbits in the Sun's Oort comet cloud), rather than 
on unbound orbits. 

In all ensembles, after approximately 10 7 yr, the mean 
number of planets per system (N p i) reaches an aver- 
a ge value between 1 . 8 and 3, similar to the findings 
of Adams fc Laughl iii ( 2003) and lPapaloizou fc Terquerril 
(|2001l ), despite the substantiall y different initial condi- 
tions (e.g., the spherical shell of iPapaloizou fc Terqueml 
l2001fh Although (N p i)(t) begins to level off at his point, 
indicating an end of strong dynamical evolution it does 
not do so entirely. A modest power law extrapolation 

6 c03s00 ensemble is an exception, with an ejection fraction of 
5%. 
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Fig. 2. — The final eccentricity distributions of the simulated 
ensembles. The colored bands show the final eccentricity distribu- 
tions of one ensemble classified as "inactive" (top), three classified 
as "partially active" (middle) and four "active" ensembles (bot- 
tom). The widths of the bands illustrate the la Poisson uncertainty 
in the distribution due to the finite sizes of simulated ensembles. 
The histograms show the observed distr ibution of eccentric ities of 
extrasolar planets with P > 20 d from Butler et al. (2006j), with 
its typical error bars shown on the left. The bin size is Ae = 0.05. 
Note the similarity of the final distributions of partially active and 
active ensembles to the observed eccentricity distribution; the ex- 
ces s of o bserved planets for e < 0.2 in the bottom panel is discussed 
in EI3.2I Ovcrplotted as a smooth solid line in the bottom panel is 
a Schwarzschild distribution with cr e = 0.3. 

near t — 10 s yr still predicts a continuing decay rate of 
~ 10% of planets per decade of time. 

Despite this slow continuing decay, after about 2 x 
10 7 yr there ceases to be any significant change in the dis- 
tributions of the orbital elements of the surviving plan- 
ets (eccentricity, inclination and semimajor axis). After 
10 s yr, for the purposes of this paper, the simulated en- 
sembles can be considered to have reached an equilibrium 
configuration. 



3. THE ECCENTRICITY DISTRIBUTION 

3.1. Classification of outcomes 

The panels of Figure [2] show the final (t — 10 8 yr) ec- 
centricity distributions of the simulated ensembles (Ta- 
ble Q]) . To visualize the expected variance due to the fi- 
nite number of systems in ensemble realizations, they are 
plotted as ±1<t wide colored bands (where Oi — \fWl/N , 
N is the total number of planets in the ensemble at 
t = 10 8 yr, and Ni is the number of planets in the i th 
bin) . The solid histogram overplotted on each of the pan- 




FlG. 3. — Initial eccentricity distributions of simulated ensem- 
bles. The distributions were obtained by binning the eccentricities 
of planets in the ensembles of Table ITl in bins of Ae = 0.05. Each 
vertex on the plot denotes the fraction of all planets in the eccen- 
tricity bin centered on that vertex. For clarity, the vertices are 
connected by straight lines. Dashed, dotted and solid lines corre- 
spond to the initial eccentricity distributions of ensembles whose 
final distributions are shown in the top, middle and bottom panels 
of Figure [2] respectively. 

els shows the distribution of eccent ricities of 125 extra - 
solar planets having P > 20 d from lButler et all (|2006f> . 
with the period condition imposed to exclude orbits that 
may have undergone tidal circularization. In Figure^ we 
plot the corresponding initial eccentricity distributions. 
The dashed, dotted and solid lines in Figure [3] denote en- 
sembles whose final distributions are shown in the top, 
middle and bottom plots of Figure respectively. 

The division into panels in Figure [2] reflects the classifi- 
cation of the ensembles as either inactive (top), partially 
active (middle), or active (bottom). This subjective di- 
vision is based on two qualitative criteria: (i) the mu- 
tual likeness of the final eccentricity distributions (easily 
seen in case of the active ensembles), and (ii) the degree 
of evolution away from the initial conditions (character- 
ized by the fraction of planets per system surviving to 
t = 10 8 yr, or the change in shape from the initial to 
the final eccentricity distribution) . The hope is that this 
classification reflects the degree of dynamical relaxation 
that the ensembles have experienced. 

The first criterion clearly separates the four "active" 
ensembles from the others (Figure[2l bottom panel). The 
eccentricity distributions of these ensembles look very 
similar, going to zero at e = and e = 1 and peaking 
around e ~ 0.3. Quantitatively, the medians (0.35 < e < 
0.4) as well as the widths 7 (0.15 < SIQR < 0.16) of the 
eccentricity distributions are all within a narrow range of 
values. The agreement is particularly striking for ensem- 
bles cl0u80 and c50s05, considering their substantially 
different initial conditions (Figure [3]) — note in partic- 
ular that in c50s05 the typical eccentricities grow, while 
in cl0u80 they shrink, yet both converge to a common 
distribution. The classification of these ensembles as ac- 
tive is also consistent with the second criterion, as they 
have all undergone substantial dynamical evolution (e.g., 
as evidenced by the reduction in mean number of planets 
per system, Figure [T] and Table [2|) . All of this indicates 
that the final eccentricity distribution of these systems 
does not retain much memory of the initial conditions, 

7 As measured from the semi-interquartile range, SIQR. 
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and is primarily a result of dynamical relaxation. Its 
general features are well described by the Schwarzschild 
distribution of equation (J2J) with a e = 0.3, which is over- 
plotted as a solid curve in the bottom panel of Figure [2j 

The second criterion sets apart the c03s00 ensemble, 
the only one classified as "inactive". Its initial and fi- 
nal eccentricity distributions are very similar, strongly 
peaked at e = 0, and inspection of the individual sys- 
tems in the ensemble shows that most have not changed 
significantly from their initial state by the end of the in- 
tegrations at 10 s yr. In the context of this paper, such 
ensembles are uninteresting as Stage 2 evolution has little 
effect on their properties. 

We assign the remaining three ensembles to an inter- 
mediate "partially active" class. These ensembles have 
undergone substantial evolution, but their final eccentric- 
ity distributions are unlike those of the active systems. 
The difference is most pronounced at low eccentricities 
(e < 0.3), while the high-eccentricity (e > 0.3) tails are 
similar to those of the active ensembles. The described 
behavior could be understood as a consequence of partial 
relaxation. For example, it could be that there is a range 
of relaxation times in different systems in a partially ac- 
tive ensemble, so that some of the systems relax fully, 
while the rest retain some memory of their initial con- 
ditions. Another possibility is that the high-eccentricity 
tail of the distribution settles into equilibrium rapidly, 
while the equilibrium form for e < 0.3 is established more 
slowly. The partially active distributions can also be re- 
garded approximately as a superposition of the distribu- 
tion of fully relaxed systems, ~ S(e; 0.3), and the distri- 
bution specified by the initial conditions. For our partic- 
ular choice of the initial conditions, the former dominates 
the high eccentricities (e > 0.3), while the latter domi- 
nates for e < 0.3. 

In Sj5j we will return to the question of this classifica- 
tion and attempt to justify it further in a more quanti- 
tative manner. 

3.2. Comparison to observations 

We compare the simulations with the distributi on of 
eccentricity of all planets in the Butler et alJ (|2006t ) cat- 
alog with period P > 20 d. This catalog reflects strong 
selection effects, primarily in mass and semimajor axis. 
Selection effects should not strongly affect the eccentric- 
ity distribution since the eccentricities of the simulated 
planets correlate with neither mass nor semimajor axis 
(this will be shown in § §3.51 and 13.6(1 . and since the se- 
lection effects in eccentricity are relatively small, at least 
in a region (P < 5yr, e < 0-6) tha t contains most of the 
observed planets (jCum ming 2004]). 

Partial qualitative agreement exists between the eccen- 
tricity distributions of the four active and three partially 
active ensembles and the observed extrasolar planet dis- 
tribution, shown in the bottom panel of Figure [2 De- 
spite this approximate agreement, there are obvious dif- 
ferences between the observed eccentricity distribution 
and the distribution of active systems: (i) an excess of ob- 
served low-eccentricity planets, and (ii) a (small) deficit 
of observed planets with high eccentricities. The two 
are inter-related, since the distribution functions /(e) are 

normalized so that J Q /(e) de = 1. 
The excess of low-eccentricity observed planets for 
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Fig. 4. — The cumulative distribution functions of eccentric- 
ity in the observed sample (solid line with no symbol) and the four 
active ensembles (solid lines with symbols), and the cumulative dis- 
tribution function of the Schwarzschild distribution S(e; 0.3) (eq.[2] 
dotted line) for all planets (top) and planets with e > 0.2 (bottom). 

e < 0.2 in the bottom panel of Figure [2] arises because 
in this region the observed systems have /(e) ~ const 
while the simulated systems have /(e) oc e. Note that 
the linear dependence seen in the simulated systems is 
a generic consequence of dynamical relaxation: in a re- 
laxed distribution, there is no reason to expect a singu- 
larity in phase-space density at e = and therefore the 
dependence /(e) tx e is simply due to the smaller vol- 
ume of phase space available near the origin (the radial 
action J r oc e 2 at fixed semimajor axis so dJ r oc ede). 
Thus the excess of observed planets at low eccentricities 
is presumably due to systems that have not experienced 
the period of dynamical activity that our active-ensemble 
simulations describe. 

On the other hand, the simulations are highly success- 
ful in reproducing the mid- and high-eccentricity part 
of the observed distribution (e > 0.2). They all success- 
fully reproduce the peak of the observed distribution near 
e ~ 0.3, as well as the decline towards e = 1 and the gen- 
eral shape of this decline. They appear to predict some- 
what more high-eccentricity planets than are observed, 
even when the e < 0.2 excess is taken into account, but 
the excess is only marginally statistically significant and 
also is consist ent with the effec ts of observational bias. 
In particular, ICummingl (|2004D showed that while the 
detection efficiency of radial velocity surveys is roughly 
constant for planets with e < 0.6, it drops rapidly be- 
yond e ~ 0.6. We have conducted simple tests, in which 
we apply a linearly or quadratically decreasing detection 
efficiency 8 to simulated data at e > 0.6, and find that 
the correction works in the direction of improving the 
agreement with our simulations. But overall, its effect 
on the shape of the observed distribution is negligible for 
the current data. 

We quantify the discussed similarities and differences 

8 Pdet = 1 — x and Pdet = 1 — x 2 , where x = (e — 0.6)/0.4. 
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TABLE 3 

KS TESTS OF ECCENTRICITY DISTRIBUTIONS (ACTIVE ENSEMBLES) 



All 


cl0u80 


cl0s40 


c50s05 


Observed 


S(e; 0.3) 


cl0s30 


0.99 


0.42 


9 x 10~ 4 


5 X 10 -3 


0.41 


cl0u80 




0.21 


3 x 10 — 5 


2 X 10 — 3 


0.15 


cl0s40 






3 x 10~ 3 


5 x 10~ 4 


8 X 10~ 4 


c50s05 








5 x 10" 7 


2 X 10" 5 


Observed 










2 x 10" 4 


e > 0.2 


clOuSO 


cl0s40 


c50s05 


Observed 


S(e;0.3) 


cl0s30 


0.99 


0.57 


0.28 


0.54 


0.21 


cl0u80 




0.57 


0.14 


0.30 


0.02 


cl0s40 






0.48 


0.13 


1 x 10~ 4 


c50s05 








0.07 


2 x KT 4 


Observed 










0.53 



Note. — The p-values of Kolmogorov-Smirnov tests between 
the eccentricity distributions of active ensembles (cl0s30, cl0u80, 
cl0s40 and c50s05), the observed sample of extrasolar planets 
with P > 20 d ("Observed"), and the Schwarzschild distribution 
with <r e = 0.3 (last column). When testing against the S(e;0.3) 
distribution the null hypothesis is that the sample was drawn 
from it (one-sample KS test). Everywhere else, the null hypoth- 
esis is that the pair of samples being tested was drawn from the 
same (but unspecified) underlying distribution (two-sample KS 
test). The top set of rows shows the results when planets of all 
eccentricities are included in the test. The bottom set shows the 
results for the subsample with e > 0.2. 



TABLE 4 

KS TESTS OF ECCENTRICITY DISTRIBUTIONS (PARTIALLY 
ACTIVE ENSEMBLES) 



All 


clOslO 


clOsOO 


Observed 


S(e;0.3) 


nlOslO 


0.11 





0.15 





clOslO 




7 x lO" 6 


0.32 





clOsOO 






0.01 





Observed 










e > 0.2 


clOslO 


clOsOO 


Observed 


S(e;0.3) 


nlOslO 


0.75 


0.55 


0.72 


1 X 10" 3 


clOslO 




0.55 


0.80 


0.06 


clOsOO 






0.68 


0.01 


Observed 








0.53 



Note. — Analogous to Table|3] but for partially active 
ensembles. 



using Figure H] and Tables [3] and [4j Figure H] com- 
pares the final cumulative eccentricity distributions of 
the active ensembles, the observed distribution, and the 
Schwarzschild distribution S(e; a e = 0.3) which captures 
the general features of the simulated distributions. The 
comparison is made both over the entire range of eccen- 
tricities (top panel), and restricted to the range e > 0.2 
to exclude possible contamination by inactive systems in 
the observed sample (bottom panel). Table [3] shows the 
p-values of the Kolmogorov-Smirnov 9 statistic calculated 
between pairs of all active ensembles, the observed dis- 
tribution, and the S(e; a e = 0.3) distribution. Table |4] 
shows the same statistic for partially active ensembles. 
At a 5% significance level 10 the eccentricity distri- 

9 The same tests were repeated using Pearson's y 2 statistic and 
led to qualitatively similar conclusions. However, for these data 
sets the x 2 statistic has less distinguishing power than the KS 
statistic (typical p-values were a factor of two larger). 

10 In this paper we adopt the a = 0.05 significance level as the 
threshold for rejecting the null hypothesis that two samples were 
drawn from the same distribution. That is, if at least 5% of pairs of 
samples randomly drawn from a given distribution would differ by 
more than the observed amount as measured by their KS statistic, 



butions of three of the four active ensembles (cl0s30, 
cl0u80, cl0s40) are pairwise 11 consistent with being 
drawn from the same eccentricity distribution (Table [3]), 
despite their quite different initial conditions (Table [1). 
The smallest p-values are obtained in tests involving 
c50s05, which began with 50 planets per system, com- 
pared to 10 for the other active ensembles. This dif- 
ference may point to a dependence of the finer details 
of the final eccentricity distribution on the initial condi- 
tions. When restricted to the subsample with e > 0.2, 
the distributions of all active ensembles are pairwise con- 
sistent. Analogous tests of the partially active ensembles 
(Table 0]) reveal similar results: substantial differences 
in the eccentricity distributions of the full samples but 
consistency among the subsamples with e > 0.2. 

Using the same statistic, we compare the final e > 0.2 
distributions of active and partially active ensembles 
with the observations (Tables [3] and [U the columns and 
rows labeled "Observed"). In all cases, the simulated 
e > 0.2 distributions and the observed distribution are 
consistent with being drawn from the same underlying 
distribution at the 5% significance level. We interpret 
this as evidence that the high-eccentricity component 
(e > 0.2) of active and partially active ensembles is pop- 
ulated by planets from dynamically relaxed systems, as 
is the majority of the high-eccentricity component of the 
observed extrasolar planet population 12 . Taking into ac- 
count the simplifying approximations of our simulation 
(giant planets only, no debris, gas disk or any other influ- 
ences, no binary companions), the exclusion of all non- 
gravitational effects (e.g., tidal effects), and the likely 
differences between our assumed initial distribution of 
masses and orbital elements and the actual distribution, 
the agreement obtained for e > 0.2 is quite remarkable. 

The Schwarzschild distribution S(e;er e = 0.3) provides 
an approximate qualitative representation of the eccen- 
tricity distribution of the active ensembles (Figure[2j bot- 
tom panel) but is not quantitatively consistent with sev- 
eral of them according to the KS statistic — probably we 
could improve the consistency by fitting er e to the ec- 
centricity distributions, but this would provide only an 
illusion of greater accuracy. 

In the scheme introduced in £13.11 the ensemble of ob- 
served planets would be classified as partially active. Its 
distribution of eccentricities may be decomposed into two 
components. One, resulting from dynamical relaxation, 
dominates the e > 0.2 regime, contains 75% of P > 20 d 
(or 55% of all) planets, and agrees well with the e > 0.2 
distributions of the active and partially active ensembles 
in our simulations. The other contains 25% of P > 20 d 
(or 45% of all) planets and dominates the e < 0.2 regime; 
here the eccentricities were set by other processes, possi- 
bly the (unknown) initial conditions or damping by low- 
mass planets, planetesimals, or residual gas. With this 

we conclude that the hypothesis cannot be rejected. 

11 We use the term "pairwise" to emphasize that these are two- 
sample tests, and that they do not test whether all ensembles are 
simultaneously consistent with being drawn from the same under- 
lying distribution. For example, three ensembles A, B and C with 
two-sample KS statistic p-values p(AB), p(BC), and p(AC) above 
a significance level a may produce p(ABC) < a in a three-sample 
test. 

12 Othe r mechanisms, such as Kozai oscillations due to a 
companion [|Wu fc Mur ray 2003; Fabrv ckv fc Trem ainc 200'3), may 
be responsible for a minority of high-eccentricity planets. 
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Fig. 5. — Time evolution of eccentricity distributions of active 
ensembles. The top, middle, and bottom panel show the eccen- 
tricity distributions at t = 10 5 , 10 6 and 10 7 yr respectively. The 
meaning of symbols is the same as in Figure [2] where the corre- 
sponding distributions at 10 s yr are shown in the bottom panel. 

decomposition in mind, for the remainder of this paper 
we will mostly concentrate on the properties of the re- 
laxed component (the active ensembles). 

3.3. Time evolution 

The panels of Figure [5] show the time evolution of 
the eccentricity distributions of active ensembles from 
t = 10 5 yr (top) to t = 10 7 yr (bottom). At t = 10 5 yr 
the distributions are still dissimilar, especially that of the 
c50s05 ensemble which is still undergoing strong dynami- 
cal activity (a consequence of the larger initial number of 
planets). By t = 10 6 yr the fraction of high eccentricity 
planets has been reduced in all ensembles, with a simul- 
taneous increase in the frequency of planets of moderate 
eccentricity. At t — 10 7 yr the eccentricity distributions 
have largely converged to a common characteristic shape, 
with the biggest change from 10 7 to 10 s yr being a fur- 
ther reduction in the number of high eccentricity planets, 
primarily by ejections (see bottom panel of Figure [2]). 

3.4. Influence of collisions 

By comparing ensembles nlOslO and clOslO we tested 
the influence of planet-planet collisions on the shape of 
the final eccentricity distribution. These two ensem- 
bles share the same initial conditions, except that in the 
case of clOslO collisions were allowed to occur, while for 
nlOslO they were not (i.e., the planets were assumed to 
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Fig. 6. — Comparison of eccentricity distributions of planets 
having a < 1 AU at t = 10 s yr, in ensembles nlOslO and clOslO. 
The meaning of symbols is the same as in Figure [2] 

have density p — 1 gcm~ 3 in the first case, and to be 
point masses in the second). We found no significant 
difference in the outcomes of these two cases (Figure [3J 
middle panel), and the final eccentricity distributions are 
consistent with being drawn from the same distribution 
at the 10% significance level (Table 2]). Since the effects 
of planet-planet collisions are likely to be most noticeable 
at small semimajor axes, we also compared the distribu- 
tions of nlOslO and clOslO planets having semimajor axis 
a < 1 AU (Figure [5]) , but again found no significant dif- 
ference. 

This outcome was not unexpected, since we have al- 
ready observed that the final eccentricity distribution is 
established over long timescales (~ 10 5 -10 8 yr, Figure[5]), 
while planet-planet collisions are infrequent events (Ta- 
ble[2]) which preferentially occur early (t < 10 4 yr). Phys- 
ically, the ratio of the escape speed from the planets to 
the typical encounter speeds between planets at 1 AU is 
large enough that the planets act like point masses. 

3.5. Dependence on semimajor axis 

Figures [7] and [5] examine the semimajor axis distribu- 
tion and the correlation of eccentricity with semimajor 
axis in the four active ensembles. 

The initial semimajor axis distribution of all of our en- 
sembles was uniform in log a. The final distribution in 
the active ensembles shows depletion at low a and some 
spreading beyond 100 AU (the initial outer limit). Nei- 
ther is particularly surprising given that interactions are 
strongest at small semimajor axes (thus the depletion at 
low a) , and that planet-planet scattering tends to spread 
out the semimajor axis distribution. The efficiency of 
depletion is particularly striking in case of c50s05, where 
fewer than 2% of planets remain on orbits closer than 
a = 1 AU; nevertheless, in general it appears that the 
semimajor axis distribution in active systems retains a 
strong memory of the initial conditions. 

Of greater consequence for comparisons of the eccen- 
tricity distributions with the observations is the cor- 
relation of eccentricity and semimajor axis. We find 
no significant e-a correlation, except for a small subset 
(less than 5% of the total) of planets scattered to high- 
eccentricity orbits at a > 100 AU (Figure [5]); this corre- 
lation is expected, since such orbits are likely to result 
from close encounters with planets at smaller radii, and 
hence must have pericenters q < 100 AU so a = q/(l — e) 
is correlated with e. The widths of the eccentricity distri- 
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Fig. 7. — The distribution of semimajor axes in the four active ensembles. The semimajor axes of known exoplanets having P > 20 d 
are plotted as a histogram. Note that the two distributions are not directly comparable, due to the significant observational biases existing 
in the latter. The bin size is A log a = 0.2. 



butions at fixed a remain approximately constant (semi- 
interquartile range SIQR~ 0.15) over the entire range 
of semimajor axes, indicating that the shape of the ec- 
centricity distribution does not appreciably vary with a 
either. This has already been implicitly assumed in £|3 . 21 
where we compared the distributions of the entire sim- 
ulated sample (with median a ~ 7.5AU-8.5AU in clO- 
ensembles and a ~ 34 AU in c50s05) to the observed 
sample (a ~ 1.3 AU). 

3.6. Mass-eccentricity correlations 

In Figure [5] we show the final (t = 10 8 yr) distri- 
butions of Ms'ml in the four active ensembles (solid 
lines) , compared to the observed M sin / distribution 
(solid histograms) and the initial mass distribution func- 
tion (dN oc M^dM, dotted line). To obtain MsinJ 
from simulations, we assume that the orbit normals are 
uniformly and randomly distributed on a sphere and as- 
sign the inclinations accordingly. Note that the observed 
distribution of exoplanet M sin / is heavily biased by the 



difficulty of detecting low-mass planets; the true distri- 
bution almost certainly is steeper and extends to lower 
values than the measured one. 

The mass distributions for the three ensembles that 
began with N p i = 10 planets per system (cl0s40, cl0s30, 
and cl0u80) converge to a similar final shape by t — 
10 8 yr. The Msini distribution of c50s05 (started with 
N p i — 50 planets per system) is different: strongly 
peaked around MsinJ ~ 7.5Mj, with an extended tail 
beyond MsmI = lOMj due to a significantly higher 
fraction of mergers than the other ensembles (Table [5]) . 

Relative to the initial conditions (dotted line), all of 
these simulations show a sharp reduction in the fraction 
of planets with small masses (M < Mj), arising because 
low-mass planets are more readily removed from the sys- 
tems, and mergers shift the distribution towards higher 
masses. Compared to the observed distribution, these 
three ensembles show a statistically significant deficit of 
planets in the range Mj < MsinJ < 3Mj but agree 
with the observed distribution for MsmI > 3Mj (KS 
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Fig. 8. — The correlation of semimajor axis and eccentricity of 
active ensembles (top four panels) and the observed (P > 20 d) 
sample (bottom panel). The error bars show the semi-intcrquartilc 
range. 

test; 5% significance level). The difference at low masses 
is not due to selection effects in the observed sample, 
since these preferentially delete the lower mass planets, 
not the higher mass ones. Simulations using a steeper 
initial mass distribution may result in a better fit to the 
observational data. 

Figure \W\ shows the correlation of eccentricity and 
Ms'ml at t — 10 8 yr in the four active ensembles and in 
the observed sample. We find non-zero but statistically 
insignificant positive mass-eccentricity correlations in the 
simulated ensembles. The median eccentricity of planets 
in different mass bins is roughly constant (e~ 0.35), as 
is the width esiQR ~ 0.15 of the eccentricity distribu- 
tions as measured by the semi-interquartile range. The 
lack of strong correlation between e and M sin I may 
be surprising given that one might expect some kind of 
equipartition, in which the most massive planets acquire 
the smallest eccentricities in planet-planet scattering. In- 
deed, as we show in Figure [Til such correlations are 
present at some times during the simulation, although 
they are never as strong as equipartition of radial ener- 
gies would predict. 

In the panels of Figure [TTJ we show the evolution of 
the mass-eccentricity correlation for one ensemble (the 
other active ensembles exhibit similar behavior). The 
top left panel shows the dependence of average eccen- 



tricity on mass at t = 0, and the subsequent panels (left 
to right, top to bottom) show the e vs. M sin / correlation 
at t = 10 5 , 10 6 , 10 7 , and 10 s yr. The correlation of mass 
and eccentricity changes during the integration. Initially 
(0 < t < 10 5 yr), the median eccentricity of low-mass 
planets grows, with the median eccentricity of high-mass 
planets also growing but by a smaller amount. This is 
followed by a period of decline of the median eccentricity 
of both low- and high-mass planets (10 5 < t < 10 7 yr 
in Figure 111)) , with eccentricities evolving to a mass- 
independent median value e ~ 0.35 at t > 10 7 yr. 

This behavior is a consequence of dynamical interac- 
tions in the system. At early times (0 < t < 10 5 yr), 
the low-mass planets are easily excited to higher eccen- 
tricity orbits by their massive counterparts. As these 
high-eccentricity planets are gradually removed from the 
system (through close encounters, ejections, or collisions 
with the star) the average eccentricity of the remaining 
planets decreases. This is supported by the finding that 
~80% of planets that are excited to e > 0.6 are removed 
from the system by 10 s yr. The median eccentricity of 
the high-mass planets initially increases, and then grad- 
ually decreases over time through the removal of planets 
on high-eccentricity orbits and the damping of eccentric- 
ity by lower mass planets. These processes continue un- 
til enough planets are removed from the system and the 
orbits of the remaining planets become sufficiently sepa- 
rated, thus rendering the system stable for the remainder 
of the integration. 

The observed sample (Figure ITTTj bottom panel) shows 
a positive correlation between mass and eccentricity at 
the ~ 5.5a level. This signal comes largely from a dif- 
ference in median eccentricity between the planets with 
Afsin/ < Mj and those with Msin/ > 3Mj. If the 
sample is restricted to the subset with M sin/ > Mj, 
its significance drops to ~ 2.8a . Some of this correla- 
tion may result from selection effects (planets are harder 
to detect either if the mass is low or the eccentricity is 
high), so we prefer to wait for more data before inves- 
tigating the implications of a possible mass-eccentricity 
correlation. 

4. INCLINATIONS OF ACTIVE SYSTEMS 

Cumulative and differential distributions of inclina- 
tions in active ensembles are shown in Figure 1121 The 
left column shows the distribution of inclinations rela- 
tive to the symmetry plane of the initial conditions (Ir), 
while on the right the inclinations are computed relative 
to the invariable plane (the plane perpendicular to the 
total angular momentum) of each system at the end of 
the simulation (7). 

The inclinations I referred to the invariable plane are 
in principle measurable by precision astrometry of multi- 
planet systems. Such measurements are currently out 
of reach of ground-based telescopes, but should become 
feasible with the launch of the SIM PlanetQuest mis- 
sion 13 . On the other hand, if the symmetry axis of the 
initial conditions is identified with the axis of stellar ro- 
tation, the inclinations Ir can be identified with stellar 
spin-planetary orbit misalignments A and are measurable 
for transiting planets usin g the Rossiter-McLaughlin ef- 
fect (the RM effect; e.g.. lOhta etal] 120051 : IWinn et all 

13 http:/ /planetquest.jpl. nasa.gov/SIM/sim_index.cfm 
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Fig. 9. — Distribution of M sin i" at t = 10 8 yr in the four active ensembles (solid line), compared to the observed distribution of M sin / 
(solid histogram), and the distribution of initial conditions (dN oc dlogM, —1 < log(M/Mj) < 1; dotted line). To obtain MsinI in the 
simulations we assume that the orbit normals are uniformly distributed on a sphere. The bin size is Alog(MsinZ) = 0.1. 

semimajor axis (Figure [T5]). in the sense that the inclina- 
tions of inner planets are on average higher than those of 
the outer. The strongest dependence is seen for ensemble 
cl0s30 (dlu/dloga — —2.8 deg/dex), while the effect is 
weakest for cl0s40 {dl^/dloga = —0.9 deg/dex). 

Until recently, all measurements of the projected an- 
gle A between the stellar spin axis and the plane- 
tary orbit axis from the RM effec t were either small 
(5, 5°) or cons i stent with zero rtQueloz et al.l 1200... 
Winn et al.ll2TM [20061 : iNarita et al.ll2007afc iWolf et al 



2005). Other mechanisms, such as Kozai oscillations 
plus tidal friction, can also ca use spin-orbit misalignment 
i|Fabrvckv fc Tremainell2007f) . 

The final inclination distributions of the three active 
ensembles that start with 10 planets appear similar, 
at least for I R > 4°, with medians 7° < I R < 9° 
(4° < 7 < 6°). For the 50 planet ensemble c50s05, the in- 
clinations are larger (median Ir = 19° and / = 10°) and 
the shape of the distribution is different. All ensembles 
have a significant fraction of planets at high inclinations 
at the end of the integrations; 10% of planets of the clO- 
ensembles possess inclinations Ir > 25° (/ > 20°), while 
10% of c50s05 planets are inclined by I R . > 51° (I > 40°). 
There is no strong correlation of inclination and eccen- 
tricity (Figure [T3|) except for the most eccentric planets 
(e > 0.7). A weak correlation exists with mass (Fig- 
ure Q3J), in the sense that the inclinations of less massive 
planets are more easily excited than those of more mas- 
sive ones. Inclinations are weakly correlated with the 



2007T ) . While most of these involved (possibly tidally 



evolved) hot Jupiters, no misalignment was found even 
in t he case of the significantly eccentric HAT-P-2b (e = 
0.5; iBakos et al.l 120071: IWinn et all 120071 : lLoeillet et all 
120071 ). A recent exception is HP 17156b (P = 21 .2 d, 
e = .67. iFischer et al.l [20071 iBarbieri et al.l l2007f ). for 
which INarita et alJ (|2007hl ) report a possible spin-orbit 
misalignment of A = 62° ± 25°. At the nominal value 
A = 62° this planet would be unusual in the context 
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Fig. 10. — Dependence of eccentricity on Afsin/ at 10 8 yr 
for the four active ensembles (top four panels) and the observed 
(P > 20 d) planets. Dots mark the eccentricity and MsinI of 
each planet. To obtain Msin/ in the simulations we assume that 
the orbit normals are uniformly distributed on a sphere. The me- 
dian and semi-interquartile ranges of e, calculated in bins of size 
Alog(Af sin/) = 0.25, are marked by the solid squares and error 
bars, respectively. The values of the Spearman rank correlation 
coefficient r s evaluated over all planets in the sample and the cor- 
responding normally distributed N (fi = 0, a = 1) variable t B are 
given in the upper right corner of each panel. The number in 
parentheses in the bottom panel is the value of t s for a subsample 
satisfying M sin / > Mj. 

of our simulations (a 1/100 event in the clO- ensembles 
and a 1/20 event in the c50s05 ensemble), even when 
the inclination-scmimajor axis correlation (Figure [15]) is 
taken into account. However, if A is lower by la then 
the inclination of HD 17156b lies in more plausible range 
(p(I R > 37°) > 0.05 for the clO- ensembles, and - 0.2 
for the c50s05 ensemble). 

Collecting a larger sample of accurate measurements 
of the projected spin-orbit angle A for eccentric planets 
may prove to be a useful endeavor. The measurement of 
a significant misalignment in HD 17156b suggests that 
misalignments are common; more comprehensive simu- 
lations of Stage 2 evolution can produce firm predictions 
for the dependence of the final inclination distribution 
on the initial conditions (initial number of planets, mass, 
semimajor axis, and inclination distribution, etc.); and 
other processes such as Kozai oscillations yield equally 
firm but quite different predictions. 



5. A MEASURE OF DYNAMICAL ACTIVITY 

The eccentricities of the three active ensembles with 10 
initial planets, as well as the e > 0.2 subsamples of all ac- 
tive ensembles, are pairwise consistent with being drawn 
from the same underlying distribution. The same holds 
true for the e > 0.2 subsamples of all partially active 
ensembles. In $3] we have taken this agreement as evi- 
dence that these ensembles have converged to the same 
eccentricity distribution. We hypothesize that this dis- 
tribution, empirically described by a Schwarzschild dis- 
tribution with er e ~ 0.3 (eq. [2]), is the equilibrium end- 
point of "dynamically active" planetary systems, where 
by "dynamically active" we mean systems whose planets 
experience strong mutual interactions and frequent en- 
counters. We now attempt to find an empirical measure 
of whether a planetary system with given initial condi- 
tions will be dynamically active. 

In the restricted three-body problem, the natural mea- 
sure of the radius of influence of a planet on a nearby test 
particle in a nearly circular orbit is its Hill radius: 



R h (r,M)=r 



M 
3Mp 



1/3 



(3) 



where M and r are the mass and orbital radius of the 
planet while M Q is the mass of the star. In the case of two 
bodies with masses Mi and M2 that are small compared 
to Mq, on nearly circular orbits with similar radii r, the 
Hill radius is ob tained from equation ([3]) by replacing M 
with Mi + M 2 (jHenon fc PetitJll986h . For the purposes 
of this paper, where we must deal with planets having 
different masses and orbital radii, we define the mutual 
Hill radius 



RH = \[Rh{r A ,M A ) 



R h (r B ,M B )} 



(4) 



as the average of the Hill radii of the two individual plan- 
ets. This definition is somewhat arbitrary but reduces to 
the usual Hill radius when one planet is much more mas- 
sive than the other. 

In the case of th e g eneral three-bod y problem 
iMarchal fc Bozisl (|1982D and lGladmar] (|1993h have shown 
that two small planets on circular, nearly coplanar or- 
bits can have no close encounters (are "Hill stable") 
if their semimajor axes are separated by 14 ai — a\ > 
2\f3Rh[\{ai + 02), Mi + Ma]. For systems of three or 
more planets, this criterion still approximately predicts 
whether the system is unstable on short timescales. How- 
ever, it does not accu rately predict long-term stability 
(|Chambers et al.lll996l) . Nevertheless, since it is usually 
the case that only the one or two closest neighbors (ex- 
pressed in Hill radii) are responsible for most of the evo- 
lution, we can use the concept of Hill radii to explore 
approximate criteria for whether a given planetary sys- 
tem is dynamically active. 

For a given planet A, we introduce the notion of its 
"Hill neighbor" B, which is the planet of larger mass 



14 Note that Gladman ( 1993) defines the mutual Hill radius as 

1/3 



Rh.g 



ai + a 2 1 mi 



■ m 2 



MI, 



O 



(5) 



In the limit of equal-mass planets on nearby circular orbits, defini- 
tions (J4j and J5} differ by a factor of \/2. 
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Fig. 11. — The evolution of the mass-eccentricity correlation in the cl0s30 ensemble. The panels (left to right, top to bottom) show the 
masses and eccentricities of planets at times t = 0, 10 5 , 10 6 , 10 7 , and 10 8 yr. The meaning of the symbols is the same as in Figure ITDI 



whose orbit comes closest, in terms of mutual Hill radii, 
to the orbit of planet A. We define the Hill neighbor 
separation, Dh, to be the minimum distance between 
the orbital ellipses of planet A and its Hill neighbor B, 
expressed in mutual Hill radii Rh- We use the minimum 
distance, instead of, say, the average distance, because 
the mutual interaction of the two planets is strongest at 
the point of closest approach. 

For example, in a system with N planets of different 
masses, there are N — 1 Hill neighbors and N — 1 Hill 
neighbor separations. In the nlOslO ensemble of Table[TJ 
there are 1000 systems of 10 planets each (at t = 0), and 
thus 9000 Hill neighbors and 9000 Hill neighbor separa- 
tions. We can observe the time evolution of this "Hill 
neighbor separation distribution" (NSD) and its statisti- 
cal properties. 

In Figure [TBI we show the NSD of ensemble cl0s40 at 
t — (top left panel) and its evolution from t = 10 4 yr 
(top right panel) through t = 10 s yr (bottom right 
panel). At t = the systems of this ensemble are very 



tightly packed, and will be unstable by virtually any cri- 
terion based on Hill radii. The ensemble reacts to this 
strongly unstable situation by removing planets through 
collisions and ejections (see Figure [lj and by redistribut- 
ing planets so as to increase the spacing between their 
orbits. As a result, the number of planets with small Rh 
decreases rapidly, and both the peak and the median of 
the NSD shift towards higher values of Rh- At t — 10 8 yr 
all neighbors are separated by more that ARh ■ 

We repeat a similar analysis for all ensembles of Ta- 
ble [JJ Figure IT71 shows the initial (left column) and final 
(right column) NSD for each ensemble, starting with the 
one "inactive" ensemble in the top row, followed by three 
we classified as "partially active" , and then the four "ac- 
tive" ensembles. A common property of all active en- 
sembles is the strong peak at small values of Dh in the 
initial distributions. For example, while all four active 
ensembles have initial median Hill neighbor separation 
Dh < 1, the lowest for a partially active ensemble (out 
of the three ensembles classified as such) is Dh = 4.4. 
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Fig. 12. — Cumulative (top row) and differential (bottom row) inclination distributions of active systems, calculated with respect to the 
reference plane of the initial conditions (In, left panels) and the invariable plane at the end of the simulation (/, right panels). The bin 
size of the differential distribution is AIr = AI = 5° . 
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Fig. 13. — Correlation of inclination and eccentricity in active 
systems. The inclinations are with respect to the reference plane 
of the initial conditions. 



Fig. 14. — Correlation of inclination and MsinI in active sys- 
tems. The inclinations arc with respect to the reference plane of 
the initial conditions. 



The final NSDs also share a number of common char- 
acteristics. All exhibit a sharp reduction in the number 
of objects at small values of Dh- This gap near Dh = 
is more pronounced in active ensembles. All examined 



NSDs peak at Dh — 12, with the NSDs of active en- 
sembles (bottom four panels) having a similar unimodal 
distribution with a strong peak at Dh — 12, a width 
ADh ~ 8 (FWHM), and an extended tail reaching to 
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Fig. 15. — Correlation of inclination and semimajor axis in 
active systems. The inclinations are with respect to the reference 
plane of the initial conditions. 
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Fig. 16. — Evolution of the Hill neighbor separation distribution 
(NSD). The panels show (left to right, top to bottom) the evolution 
of the NSD of ensemble cl0s40 at six different times. Plotted are 
the number of planets in the ensemble in bins of width ADg = 1.5, 
where Djj is the Hill neighbor separation 




much larger values of Dh- 



Fig. 17. — Comparison of initial and final Hill neighbor separa- 
tion distributions (NSDs) of simulated ensembles. Each row shows 
the initial (t = 0, left panel) and final (t = 10 8 yr, right panel) 
NSD of the ensemble. Plotted are the number of planets in the 
ensemble in each bin of width ADg = 1.5. 

Figure [18] shows the evolution of the median Hill neigh- 
bor separation Dh for the ensembles of Table [TJ The 
division into dashed, dotted and solid ensembles, corre- 
sponding to inactive, partially active, and active ensem- 
bles, follows the convention adopted for Figure [3] and 
the separation into top, middle and bottom panels in 
Figure O The four ensembles with initial Dh < 1 are 
all active; the three having 4 < Dh < 7 are partially 
active, and the only ensemble with a large initial value 
Dh — 14 is inactive. All active and partially active en- 
sembles converge to a final median Hill neighbor separa- 
tion D H ^ 12 - 14 after 10 s yr. 

We take these results to offer hope that the initial value 
of Dh may be used as a crude measure of the dynami- 
cal activity of an ensemble, and hence as a predictor for 
the classification of the final eccentricity distribution. A 
much more thorough exploration of possible initial con- 
ditions is needed before we can tell whether this hope is 
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Fig. 18. — Evolution of median Hill neighbor separation Dh 
with time. The median is plotted against elapsed simulation time 
for each ensemble. Active, partially active, and inactive ensembles 
are plotted with solid, dotted and dashed lines, respectively. 

justified. In the meantime, our best guess is that systems 
with Dh % 10 are likely to be at least partially active. 

Other statistical measures may eventually prove to be 
more useful or reliable in characterizing the dynamical 
activity of a planetary system. When devising the one 
employed here, we were guided by the criteria that it 
should reflect the level of short-term dynamical insta- 
bility present in the system (satisfied by expressing the 
distances to the closest more massive neighbor in Hill 
radii) , and that it should be applicable both to coplanar 
orbits and to orbits with significant inclinations (requir- 
ing the relatively complex definition of Hill neighbor sep- 
aration). The statistic as defined above works reasonably 
well for systems of the type existing in our simulations 
- a few planets with a limited range of masses (two 
decades). However, care must be taken when applying it 
to (and interpreting it for) other types of systems, as we 
cannot prove that it will work everywhere equally well, 
and strongly suspect that there are possible, though per- 
haps pathological, planetary systems for which it does 
not work at all. 

6. DISCUSSION 

As described in the Introduction, the possibility that 
planet-planet interactions play a significant role in ex- 
plaining the origin of the extrasolar planet eccentricity 
distribution has been discussed since soon after the first 
extrasolar planets were discovered. Most of these discus- 
sions focused either on exploring the dynamics of sim- 
plified two or three planet systems, or on following the 
dynamical evolution of planetary systems whose initial 
conditions were inspired by planetesimal accretion the- 
ory 15 . 

In this paper, we have integrated large ensembles of 
randomly constructed planetary systems over 100 Myr 
timcspans, simulating PP « 5 x 10 12 planet-periods 

15 In particular, the initial version of our paper was submit- 
ted to the arXiv prep rint server on the same day as papers by 
IChatteriee et al.l (2007) and lFord fc R asio ( 2007), who reach many 
of the same conclusions. 



(eq. [J). The output from these simulations, and future 
simulations of this type, offers a rich resource for studies 
of Stage 2 planet evolution, and here we have focused 
on only a few aspects of this evolution, in particular the 
distributions of eccentricity, inclination, and separations. 
Initially, we classify the ensembles according to their final 
eccentricity distributions, and later show that this classi- 
fication is strongly correlated with the initial median Hill 
neighbor separation Dh, in that all ensembles classified 
as dynamically active had Dh < 1 in the initial state. 
In all dynamically active ensembles that we have exam- 
ined, we obtain the same final eccentricity distribution 
for a remarkably wide range of initial conditions. This 
distribution is described by a Schwarzschild distribution 
(eq. [2]) with a e ~ 0.3. For e > 0.2 the final eccentricity 
distribution in our simulations of active ensembles agrees 
with the observed eccentricity distribution of extrasolar 
planets remarkably well. The excess of observed systems 
with e < 0.2 may reflect either a population of planetary 
systems that are not dynamically active, or eccentricity 
damping by low-mass planets, planetesimals, or residual 
gas. In the former case, comparison with the observa- 
tions suggests that about 25% of the known extrasolar 
systems with P > 20 d are inactive, and 75% active. 

We find little or no correlation of other parameters 
(semimajor axis, planetary mass, and inclination) with 
the final eccentricity, although such correlations are 
present during periods of dynamical instability early on 
in the simulation. 

This "equilibrium" distribution of eccentricity is a 
product of dynamical relaxation of an initially unsta- 
ble system. The distribution is mostly established after 
10 7 yr and remains stable to at least 10 8 yr, where our 
simulations end. A few integrations have been carried to 
longer times, and show no evidence of further evolution. 
By 10 8 yr most of the active ensembles have only 2-3 re- 
maining planets in the three decades of semimajor axis 
that we originally populated. 

We find further that in all partially active and active 
ensembles that we examined, Dh converges to a common 
value between 12 and 14. In active ensembles, the distri- 
bution of Dh converges to a common shape as well, with 
a peak at Dh — 12 and a width ADh — 8 (full width at 
half maximum). Thus both the eccentricity distribution 
and the distribution of Hill neighbor separations in ac- 
tive ensembles appear to be a common endpoint of the 
dynamical relaxation process. 

An aggressive interpretation of the similarity of the 
observed and theoretical eccentricity distributions for 
e > 0.2 is that the high eccentricities of observed plan- 
ets have arisen as an endpoint of dynamical relaxation, 
by processes similar to those seen in the simulations of 
£|3.1i long after the initial stage of planet formation and 
dispersal of the protoplanetary gas disk were complete. 
This interpretation leads to a number of interesting con- 
clusions: 

• There exists no single "right" ensemble of initial 
conditions at the end of Stage 1. Instead, there is 
a multitude of substantially different ensembles of 
initial conditions that lead to the same or similar 
final outcomes, at least for the eccentricity distri- 
bution and the distribution of Hill neighbor sep- 
arations. An important corollary is that the de- 
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tails of initial conditions are impossible to deduce 
from the "relaxed" component of the observed ec- 
centricity distribution, except to say that they are 
likely to be in the "active" regime where dynamical 
evolution is strong enough to drive the relaxation 
process. 

In a large fraction of systems the final products 
of Stage 1 planet formation must be dynamically 
active. This is in principle possible fo r both the 
planetesimal (e.g.. iKokubo &: Idalll998f) and grav- 
itational instability (e.g., iBossl [20001 ) models. Of 
course, the separation into Stage 1 and Stage 2 is 
somewhat artificial, since the initial Stage 2 evo- 
lution occurs on timescales short compared to the 
likely duration of Stage 1. 

Planet-planet scattering in active ensembles 
changes the orbital planes of planets, broadening 
the distribution of inclinations with respect to the 
symmetry plane of the initial conditions. The same 
broadening is seen in simulati o ns of three-planet 
scattering by iChatteriee et alj (|2007D . Assuming 
that the host-star spin vector is parallel to the 
symmetry axis of the initial protoplanetary disk, 
and that no other effects re-orient the spin axis of 
the star or the invariable plane of the planets (e.g., 
Tremainc 1991), this inclination distribution is de- 
tectable, at least in principle, in transiting plan- 
ets through the Rossiter-McLaughlin effect. Given 
that the distribution of misalignments depends on 
the initial conditions more strongly than the distri- 
bution of eccentricities in active systems, its mea- 
surement in a large sample of extrasolar planets 
may yield valuable information about the endpoint 
of Stage 1 and the initial conditions for Stage 2. 

The typical final number of giant planets in active 
systems is between 2 and 3, at least over the range 
of three decades in semimajor axis that we popu- 
late in the initial states (presumably other planets 
could form well outside this semimajor axis range, 
and their interactions with our simulated planets 
would be negligible). At t — 10 8 yr, on average, 
20% of active systems remain with only 1 planet, 
75% have two or three planets, and only 5% have 
four or more planets, suggesting that most extraso- 
lar planetary systems should have 2-3 giant plan- 
ets in this semimajor axis range (about one giant 
planet per decade). Consequently, extrasolar sys- 
tems with a single detected eccentric planet are 
likely to harbor at least one more planet of com- 
parable mass. Observational data show long-term 
radial velocity trends indicative of the presence of 
another planet in ~ 50% of the curre ntly known 
exoplanet systems (|Butler et al.l l2006h . However, 
an equally interesting prediction is that such sys- 
tems are also not likely to harbor more than one 
or two additional giant planets. 

These predictions are consistent with current ob- 
servational data on the fraction of multi-planet sys- 
tems. To compare our simulations directly to the 
observations, we cull them at age 10 8 yr to keep 
only those planets that produce a radial velocity 



semi-amplitude K > 10m s _1 and have periods in 
the range 20 d< P < 2500 d (the lower limit elim- 
inates hot Jupiters, and the upper limit approxi- 
mates the longest detectable orbital period). The 
culled ensembles cl0s30, cl0s40 and cl0u80 predict 
a ratio of single- to multi-planet systems ~ 86 : 14 
(84 : 16 lowest, 89 : 11 highest), in excellent agree- 
ment with the 87 : 13 ratio seen when the observed 
sample subjected to same culling. The agreement is 
only weakly sensitive to the exact choice of thresh- 
old K or the limits imposed on P. The prediction 
is less good for the ratio of semimajor axes of the 
outer (02) and inner (ax) planet in multi-planet 
systems, where the simulations typically peak at 
4 < a 2 /ai < 8, while 80% of the observed ratios 
are in 1 < 02/01 < 4 range. The distribution of 
semimajor axis ratios is not universal across the 
three ensembles either, pointing to a dependence 
on initial conditions that warrants further investi- 
gation. Finally, the c50s05 ensemble, due to the 
efficient clearing of the zone with a < 1 AU (Fig- 
ure \Ty, predicts no observed multi-planet systems 
at all. 

• The typical final separation of planets in active 
multi-planet systems should be Dh — 12-14. The 
determination of this statistic in the currently 
known multiplanet systems is made difficult by the 
unknown inclinations, both the distribution of in- 
clinations relative to the invariable plane and the 
inclination of the invariable plane to the line of 
sight. We define a new statistic D' H , which is ob- 
tained from Dh by replacing all planet masses M 
by Msin/, and compute D' H for our ensembles by 
assuming that the normal to the invariable plane is 
distributed randomly and uniformly on a sphere 
and culling the ensembles using the same crite- 
ria on period and velocity amplitude described in 
the preceding paragraph. We obtain D' H ~ 12- 
13. To compute the analogous statistic for the 
observations, we assume that the inclinations rel- 
ative to the invariable plane are described by a 
Schwarzschild distribution (Eq. [2]) with 07 = 10°, 
the distribution of the nodal longitudes fl is uni- 
form, and we take the longitudes of periastron w 
from lButler et all (|2006f ). We find D' H ~ 8 for the 
13 known multi-planet systems. This agreement 
is probably adequate given the large statistical er- 
rors, although it is also possible that the observed 
systems can be stable at smaller values of D' H than 
our simulations because some of the planets are 
stabilized by mean-motion resonances. 

• The median separation Dh in the solar system is 
15.6 (four giant planets) or 21.2 (all eight planets). 
Our results therefore suggest that the solar system 
is inactive, which is consistent with the observation 
that the eccentricities of planets in the solar system 
are much lower than in extrasolar planetary sys- 
tems 16 . Our crude estimate in £13.21 suggests that 

16 Although the solar system appears to be inactive, in the sense 
that Djj appears to be larger than needed for long-term stability 
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at least 25% of extrasolar planetary systems are 
inactive, so in this respect the solar system is not 
unusual. 

• All active systems that we have examined eject a 
significant portion of the initial mass in planets. On 
average, the active systems we simulated ejected 
50% of the initial mass (10Mj ejected) if they 
started with 10 planets, and 80% (90Mj ejected) if 
they started with 50. Therefore, free-floating plan- 
ets should be common and have a number density 
roughly comparable to the number density of ex- 
trasolar planets. Such planets may be detected 
by future microlensing surveys, or in open clus- 
ters as planetary-mass objec ts not bound to stars 
(|Zapatero Osorio et al.ll2000D . 

The constraints and limitations of the above conclu- 
sions have to be kept in mind. The distribution of eccen- 
tricities is only a one-dimensional projection of the multi- 
dimensional distribution of orbital elements, masses and 
other properties, and it is likely that other statistics of 
this distribution do not approach universal values in ac- 
tive Stage 2 evolution (for example, the mass and semi- 
major axis distributions). It is also a priori possible that 
the equilibrium eccentricity distribution may be different 
for some ranges of initial conditions that we failed to ex- 
plore in this paper. A much broader exploration of the 
possible initial conditions for Stage 2, and of the distri- 
bution of orbital elements and masses at the end of Stage 

of the existing planets, there are almost no locations in the outer 
solar system (between Jupiter and Neptune) in which additional 



2 evolution, is needed to investigate these questions. It 
is particularly important to explore (i) a steeper initial 
mass function (more low-mass and fewer high-mass plan- 
ets) , since the final mass distributions in our simulations 
probably have too few low-mass planets (Figure [9|) ; (ii) 
a minimum-mass cutoff lower than our current value of 
0.1-Mj, since a larger population of low-mass planets may 
damp eccentricities; (iii) tidal dissipation from the host 
star, which may affect orbits with pericenters less than a 
few stellar radii; (iv) active systems that initially contain 
giant planets only beyond a few AU, to investigate what 
fraction of giant planets could acquire small semimajor 
axes through planet-planet interactions; (v) active sys- 
tems containing terrestrial planets. Although investigat- 
ing the wide range of possible initial conditions and final 
states of Stage 2 evolution is a massive task, the only ma- 
jor resource required for this investigation is processing 
time on cluster computers. 
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planets could be inserted on stable orbits (Holman 1997). 
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